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The ground state energy of the 2-D Wigner crystal is determined as a function of the thickness 
of the electron layer and the crystal structure. The method of evaluating the exchange-correlation 
energy is tested using known results for the infinitely-thin 2D system. Two methods, one based on 
the local-density approximation (LDA), and another based on the constant-density approximation 
(CDA) are established by comparing with quantum Monte-Carlo (QMC) results. The LDA and 
CDA estimates for the Wigner transition of the perfect 2D fluid are at ra = 38 and 32 respectively, 
compared with = 35 ± 5 from QMC. For thick-2D layers as found in Hetero-junction- insulated- 
gate field-effect transistors, the LDA and CDA predictions of the Wigner transition are at rs = 20.5 
and 15.5 respectively. Impurity effects are not considered here. 
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I. INTRODUCTION 



Two-dimensional (2D) electron layers exist, for exam- 
ple, at the interface between GaAs and Gai-^Al^; As, 
or at the interface of a metal oxide and a semiconduc- 
tor. Such interfaces are important in field-effect transis- 
tors and other devices. The state of this 2D electrons 
system can be fluid or, at sufficiently low density, the 
electrons condense and crystallizeP, We define the 
2D-density parameter given by AjN — Trr^^ where A 
is the area occupied by the A'^ electrons. The Wigner 
solid appears for > 35ao where oq = h^/{m*e*'^) 
is the Bohr radius. Here m*,e* are effective parameters 
for the mass and the charge of the electron, and absorb 
the dielectric constant, the band mass and other mate- 
rial properties of the system. Thus in GaAs, the effective 
atomic unit of energy is reduced from 27.21 eV to the 
milli-Volt range. These (reduced) atomic units, such that 
m* = e* = h=l will be assumed through out this paper. 
There has been many studies of 2D-clcctron liquids or 
Wigner crystals 111111110,11, especially using quantum 
Monte Carlo (QMC) simulations and other methodsHQ 
assuming that the 2-D layers are infinitely thin. How- 
ever, although the 2D electrons reside in the {x, y)-plane, 
they have a transverse density 77(2) confined to the lowest 
subband of the hetero-structure [l| . While the quasi-2D 
electron liquid has recently seen much attentiori, both 
experimentaljO, ^-nd theoretical floL diii4, 13, the 
Wigner crystal in thick 2D layers has not been followed 
up since the work of Fujiki and Geldart ^1|. Fujiki et 
al., have determined the effect of the 2D-layer thickness 
on the electrostatic energy and found that the hexagonal 
lattice is the most-stable crystal structure, as with the 
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(5-thin 2D layer (52D). However, they did not consider 
the effect of exchange and correlation which is usually 
addressed via Quantum Monte Carlo methods, or via a 
detailed_analysis of the correlated phonons in the electron 
crystal 7] . We note that recent Hartree-Fock (HF) calcu- 
lations of 62D Wigner crystals using large plane- wave ba- 
sis sets, e.g, those of Trail et al0., seem to recover a HF 
energy nearly identical to the single- Gaussian harmonic 
approximation for localized electrons. Various aspects of 
such a model have been considered in a brief but insight- 
ful paper by Nagvp^Hoj. In this study we show that the 
single-Gaussian approximation, and the local-density ap- 
proximation (LDA), can recover the QMC total energy 
with surprising accuracy. We also show that a method 
based on constructing a constant-de nsity a p proximation 
(CDA) to the inhomogeneous densitv |10ll2(]| | can be prof- 
itably used for calculating the electrostatic potentials and 
the exchange-correlation energies of these systems. 

The plan of the paper is as follows. In section^lwe in- 
troduce the Hamiltonian, the effective Coulomb interac- 
tion in quasi-2D layers, and calculate the electrostatic en- 
ergy of the lattice for several 2D-crystal structures. Here 
we use the CDA to replace Fang-Howard type densities 
in the z-directionll, 'lO'|, thus simplifying the analytical 
work. The details of lattice-sum evaluations are relegated 
to an appendix. In section ITTll we consider the 5-thin 2D 
layer and present results for the gaussian-localized model, 
we also present the exchange-correlation energy E^c cal- 
culation using the CDA and the LDA. The resulting to- 
tal energy is very close to that of QMC and recovers a 
liquid^solid Wigner transition (WT) at rs ~32 to 38, 
while the current QMC estimate is rs = 35 ± 5. In sec- 
tion we consider Gaussian- localized 2D systems with 
finite thickness, for 2D layers found in HIGFETS. That 
is, in systems where the layer thickness is also defined by 
the sheet density, as in Tan et al 11, 12J. Here we have 
no QMC results for comparison. The total energy of the 
quasi-2D Wigner crystal is compared with the total en- 
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ergy of the quasi- 2D liquid T^. Here the WT is found to 
occur at ~15 to 21 in quasi-2D layers realized in clean 
HIGFETS, 



II. THE COULOMB ENERGIES OF 2D 
LATTICES 



The Hamiltonian of our system is, in atomic units, 



bb 



(1) 



where the first term is the kinetic energy of the elec- 
trons. The three remaining terms are the electron- 
electron interaction and the interactions involving the 
uniform, static neutralizing background, indicated by the 
subscript b. This neutralizing background arises from a 
homogeneous layer of donor ions which have acquired 
a positive charge after donating their valence electrons 
in forming the 2D-electron layer. We assume that the 
electron layer is confined near the plane z — and ex- 
tends into the region of z > due to the finite width of 
the envelope function. The donor ions are modeled by 
a homogeneous layer of positive charge of areal density 
Pd = N/A, situated at z = —bd, where 6^ = is a 
positive quantity. The z-direction density is r]{z), and in 
the plane, an areal density Pe(r) with r = {x,y). The 
subband distribution 77(2) is usually modeled by a Fang- 



Howard distribution rjfh{z) = (l/26^)z^ 



(n.b.. 



b — l/b used in Ref. or various other forms, e.g, that 
of a quantum well. The form of the density is obtained by 
fitting to a self-consistent calculation of the Schrodinger 
equation for the electron motion in the z-direction. In 
our work, we do not repeat this calculation, but simply 
take the value of the parameter b, or other parameters 
needed to define the self-consistent solution for the sub- 
band. Moreover, as discussed below, such inhomogeneous 
densities can be replaced by a constant-density slab hav- 
ing an equivalent electrostatic potential, using the CDA 
discussed by Dharma-wardanaj^. The CDA method 
pU] involves replacing an inhomogeneous density r]{z) by 
a slab of constant-density fj of width a linked to 77(2) by 



rj = — — / ri{z)^ dz 
a 



(2) 



This equation has also been proposed by Gori-Giorgi et 
al'2^., in a method for calculating system-adapted corre- 
lation energies. Using Eq.|21a Fang- Howard (FH) density 
of length scale b can be replace by a homogeneous density 
of width a = (16&)/3. 

Consider two electrons in a quasi-2D layer separated 
by a distance r in the 2D plane, and located at zi and 
Z2, with a FH distribution 77(2:) in the z direction. Then 
the Coulomb interaction is of the form 



W{r)^ / dzi dz 



>2 + (zi-Z2)2|l/2 



(3) 



W{r) may be written as F{r)/r where F{r) is the form 
factor. No analytic form exists if ri{z) is the FH form, 




FIG. 1: (Color online) Profiles of the Fang-Howard density 
(solid curve) for fe = 4 and its equivalent constant density 
approximation(CDA, dashed curve). Inset: the bare Coulomb 
potential 1/r and the Coulomb potential F(r)/r modifed by 
the Fang-Howard profile. The triangles are calculated using 
the CDA. The CDA width a=21.33 for 6=4 corresponds to a 
HIGFET at = 32.496. 



while the g-space form, F{q)2Tr/q is analytically avail- 
able. For GaAs/AlAs based HIGFET-like systems, it 
takes the form. 



F{q) = [1 



(4) 
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where b' = l/b and follows the definition in Ref. 0. How- 
ever, if the FH distribution is replaced by the CDA, then 
F(r) and F(q) are given by 

W{r) = V{r)F{s), s = r/a, t = y/{l + s^) (5) 
log 



and 



F{s) 



Wiq) 
F{p) 



2s 



1-t 



V{q)F{p), p^qa, V{k) 
(2/p){(e-P-l)/p+l} 



2Ti/q 



(6) 



(7) 
(8) 



The form factors F{r), F{q) are a measure of the re- 
duction of the strength of the 2D interaction due to the 
thickness effect. These results provide equivalent analyt- 
ical formulae for the FH-density, and tend to the ideal 
2D Coulomb potential when the width a tends to zero. 
Also, for HIGFETS, it is known that the FH-parameter 
b is linked to the 2D density parameter Ts- Hence it can 
be shown [l3 that 



b = (2r2/33)i/3 
a = 166/3 

/3 = = [2/(33r,)]i/3 



(9) 
(10) 

(11) 



Hence /3, the FH parameter b in units of Ts, and also 
the ratio a/r^, i.e., (z-width)/(2D-disk radius) decrease 
as Ts^'^ with increasing r^. 
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A. Coulomb energy 



where 



In the foUowing we do not at first specify the form of 
the transverse density r]{z). In calculating the Coulomb 
energy Ecou, i-e, the electrostatic energy, we isolate the 
long-range contributions which cancel in the q = limit, 
since we are dealing with a homogeneous, neutralizing, 
static background. The total Coulomb energy is the sum. 



1 r f piq(i--r') 



(13) 



Ecou = lini[Sdd(9) + E,^{q) + 2EM] (12) 

9^0 



r]{z)T]{z') 



1 /• p pOO pOO 



rjiz') 



r 



(14) 
(15) 



Edd is the interaction energy of the ions, E^e is the inter- 
action of the electron layer and E^d is the energy due to 
interaction between the ions and the electrons. To cal- 
culate these terms, we proceed as in Fujiki and Geldart 

MM. 



Edd{q) 



T^Apd^ 



N 



q qrs^ 

We introduce the integral transformation 



1 



\^ 



(16) 



(17) 



and note that 

Eed{q) d\ J d\'pdPe{v')e'^-^'-'"> dz'r^iz') 

Jo ^ 
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dz'r;(z')e-«l^'+''<^l' 



(18) 



For -Egej we use a lattice sum technique based on the 
6 Jacobi function, Ea. H19|l . and its imaginary transform, 
Ea. (pn|l . given below: 



(z,X) 



(19) 



(20) 



We decompose the lattice into rectangular sublattices in- 
dicated with sublattice vectors pj. So, the position vec- 
tors of the electrons on nodes / and J are given by 

ri = maix + na2y, rj = (m'ai + pj)x + {n'a2 + pj)y 

where m, to', n, n' are integers, oi, 02 are lattice constants 
of sublattices. For example, in a square lattice ai — a2 
and {pj} — {(0;0)}, in a hexagonal lattice 02 = \/3ai 

and {pj} ^ |(0; 0), (^^^ |. To proceed further, we 

need to specify the form of the density. If the electrons 
are assumed to be exactly localized on the nodes of the 
crystal, then 



Pesir) = ^(5(r-ri) 



(21) 



Such exact localization of the electrons provides the 
model for the classical electrostatic energy, i.e, the 
Madelung energy. In the quantum calculation we sup- 
pose that each electron is localized around a node I of 
the lattice and the wavefunction is taken to be a gaussian 
normalized over the 2D plane. 



(r) 



(22) 



The parameter a is chosen to minimize the total energy. 
Hence the localized density is 



(23) 



The gaussian- width parameter a is of the form ajr'J'^ , 
with a taking a lower-bound value of 0.5 (see Ref [T^'). 
These two forms of the density will be studied below, and 
the Gaussian approximation will be justified by compar- 
ison with results from detailed plane-wave calculations. 
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TABLE I: The Madelung energy, Ecou per electron are given 
for diflerent values of the Fang-Howard parameter /3 = b/r^ 
for hexagonal(hex), square(sq), rectangular(rec), centered 
rectangular(cr) lattices defined by their unit vectors ai:a2. 
The Ts parameter in the corresponding HIGFET, Eq. |2l is 
also given. The energies are in units of l/vs- Thus the 
Madelung energy in Hartrees for a (5-thin hexagonal lattice 



is -1.106103/r^ 


HIGFET 


<x 


60606 


60.606 


0.06060 


(al : a2) /3 -> 





10-2 


10-1 


1 


hex{V3 : 1) 


-1.106103 


-1.052959 


-0.591433 


3.144793 


ct{V2 : 1) 


-1.104080 


-1.050937 


-0.589507 


3.145401 


sq(l:l) 


-1.100244 


-1.047103 


-0.585854 


3.146555 


rec(V2:l) 


-1.078201 


-1.025072 


-0.564890 


3.153217 


rec(V3 : 1) 


-1.042843 


-0.989733 


-0.531301 


3.163948 



B. Calculation with the 5-distribution 
Using Eq. (^7)) and (HIJ we have 

E,,{q) = /'°°-^/(y)^e*q("-i-r.)e-y|ri-r.|^24) 

f{y) = / dz dz'r]{z)r]{z')e-y'^'-''^' 
Jo Jo 

The details of the evaluation are given in the Appendix. 

We have evaluated Ecou, Eq. Elfor different lattices: 
square, rectangular, hexagonal and centered-rectangular. 
The Coulomb energy depends only on f3 = {h/rg) = 
(3a)/(16rs), r = (02/01) and {pj}. Our numerical cal- 
culations of Eqou are summarized in Table ^ Results 
for /3 = 10"^ are at unrealistically low HIGFET densi- 
ties, but are of formal interest. Results for even smaller 
values of (3 may be found in Fujiki et al0,|^. A com- 
parison with the results of Ref. shows that our results 
are in agreement when a geometrical term arising from 
the slight difference in the models is taken into account. 
(As seen from the details given in the appendix, we have 
an additive term N {2a / 2>) / r s"^ = N {i2h / 9) / r s"^ in our 
calculation while Fujiki and Geldart have A^(33&/8)/rs^. 
Agreement is obtained if we replace our term by theirs). 

It is seen that the total Coulomb energy increases as 
P = b/vg increases for all cases studied. The hexagonal 
lattice has the lowest energy for all f3. Moreover, there is 
no crossing between the different energy curves for any 
of the lattice structures. 

The dependence of the total Coulomb energy of the 
centered-rectangular lattice and rectangular lattice as a 
function of the ratio r — 01/02 for the quasi-2D sys- 
tem remains similar to the 5-thin case. Two equivalent 
minima at r = y/3 and l/y/S correspond to the hexago- 
nal structure. For the rectangular lattice, the minimum 
corresponds to r = 1, i.e., to the square structure. We 
choose the range (3—0.05 to 0.5, which corresponds to 
Ts ~ 0.5 to ~500 and fit the Madelung energy of the sta- 
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FIG. 2: (Color online) (Coulomb energy per electron in 
atomic units for a perfectly 2D system and for 2D layers in a 
HIGFET, using Eq.@ to define the thickness. 

ble hexagonal lattice (see Table IJ) to the analytic form 

i=4 

i^Cou(rs,/3) = 5Zc./5Vr, (25) 

i=0 

where co = -1.106103, ci=5. 34722, C2 -2.15257, 
C3=1.48663, and C4 = -0.430473. 

In Figure El we have plotted the Coulomb energy as a 
function of rg using Eq.® to relate the thickness to the 
Ts value. We observe that the thickness of the system 
has a significant effect on the energy. 



C. Classical calculation with the gaussian 
distribution. 

If the electron distribution at each site were a gaussian, 
the classical electrostatic energy, E^e, can be calculated 
using the same techniques as before (Appendix). 

1 f Oi. \ 

2\/7rJo \y + aj 

^e"i7nj('''"'"^^'e"*i^'»'('"'"'"^) (26) 

We use the same integral separation with E^^ and 
the Jacobi function 9 and its transformation. We may 
verify that when a tends to zero, that is to say the gaus- 
sian distribution tends to the (5-distribution, E^e reduces 
to the Madelung energy of the previous section. Also, if 
there is no effective overlap among the gaussian distri- 
butions, the distributions can be replaced by equivalent 
point charges at the lattice sites and the Coulomb energy 
should reduce to the Madelung energy. However, as al- 
ready remarked by previous authors [lMll9| . the charge is 
not perfectly contained within the Wigner-Seitz disk in 
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the 2D problem. The variations of the thickness and of 
the lattice type give results similar to the (5-thin case. We 
consider the variation of a to minimize the total energy 
within a quantum calculation, and hence do not develop 
this classical calculation any further. 

III. PERFECTLY TWO-DIMENSIONAL 
SYSTEMS 

In this section we consider a perfect, i.e., S- 
thin 2D layer within a Kohn-Sham density-functional 
approachp^ to the quantum mechanics of the problem. 
Since the i5-thin 2D system has been studied extensively, 
we use it as a reference system to examine the LDA and 
the CDA as useful tools for calculations of E^c of Wigner 
crystals. The Hohenberg-Kohn theorem asserts that the 
total energy is a functional of the one-electron density, 
and that it is a minimum for the true density distribution. 
We model the one-electron density as a sum of gaussians 
centered on each lattice site, and hence the variational 
problem reduces to a determination of the width param- 
eter a of the Gaussian that minimize the total energy. 
The total energy of the system at a given can be writ- 
ten as: 

Et ^ EHF{a,rs) + E,,{a,rs) (27) 

where EHFict,'rs) is the Hartree-Fock energy of an elec- 
tron. It will be seen that this is effectively the energy of 
an electron on a single site, and moving in the potential 
well created by the gaussian distributions on other sites. 
If the gaussians were perfectly localized, the Coulomb 
energy would not depend on a. The effect of the overlap 
can be easily included in the variational problem, with 
the energy given by < iplHli/j > / < ipltp >, and this 
has an effect for small r^. Here -0 is a Slater determi- 
nant of gaussians. For the hexagonal lattice, the overlap 
contribution from two nearest-neighbour gaussians is 

Sy(r,) = exp[-(a/2)(1.09r,)2] 

where LCQr^. is the nearest-neighbour distance. Unless 
the contrary is stated, the results reported here will in- 
clude the overlap correction. The a which minimizes 
the Hartree-Fock problem is not the same as that which 
minimized the total energy inclusive of E^c- In the next 
section we look at the problem without E^c 

A. The Hartree-Fock energy Ehf 

The Hartree-Fock energy is composed of the classical 
Madelung energy which defines a constant energy term, 
plus the quantum mechanical energy associated with the 
motion of the electron in the field of the other electrons. 
Since the electrons are strongly localized, especially for 
large r^, a Slater determinant made up of one gaussian 
function at each lattice site is commonly assumed. The 



TABLE XL Comparison of the plane-wave calcuIation llTll of 
the HP energies Eh f of the (5-thin 2D hexagonal Wigner lat- 
tice with the single-Gaussian harmonic lattice energies. -EJ^^^ 
and Ehar are energies without and with the overlap correc- 
tions. 



rs - 




20 


30 


40 


60 


-Ehf 

-Ehar 
~Ehar 


X 10 
X 10 
X 10 


0.447270 
0.447155 
0.437058 


0.311642 
0.311786 
0.308344 


0.239528 
0.239822 
0.238326 


0.164036 
0.164530 
0.164113 



total energy consists of a kinetic energy term and a poten- 
tial energy term. These two terms are equal by the virial 
theorem and hence we only need to evaluate the kinetic 
energy. Usually, Hartree-Fock energies contain a sizable 
exchange contribution. However, the localized-gaussian 
exchange energy is easily shown to be negligible, and we 
called it the Wigner-exchange energy, Exwc- 

In Table ^1 we compare our localized-gaussian (har- 
monic) calculation with the results of the extensive plane- 
wave HF calculation by Trail et al|23|- The results shown 
in Tabled show that the localized single-gaussian model 
is adequate to describe the Hartree-Fock approximation 
for this svstem|23|. 

Note that our calculation is effectively an "Einstein 
model" of oscillators, and the kinetic energy is given by, 

i;i^(a) - -y < 0/|V?|0/ >=7Va (28) 

The gaussian width which minimizes the energy may be 
fitted by the form a = 0.6263/ri-57. 

Since the exchange of electrons actually involves a de- 
localization process, we believe that the exchange inte- 
gral evaluated with fixed gaussians does not lead to a 
true evaluation of the exchange in these systems. The 
Wigner-exchange energy between two electrons of spin 
Si, Sj, is by definition 

= -V^e-"('-'-'-^)',5,„.,. (29) 

We can define a polarization parameter C, = {Ni^ — 
Ni)/N, therefore 

Ex^c{a,0^-lT.Eln,c (30) 

This Exwc, niay be safely neglected for the values of a 
occurring in this problem. 

B. The CDA exchange and correlation energies 

The correlation energy is the most difficult object to 
calculate, and QMC has been the preferred approach. 
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even though this requires a a major numerical effort. 
However, for an uniform density profile, the correlation 
energy is well known Hence, as in Eq.iEl we map 
the inhomogeneous density in the {x,y) plane, p(r), to a 
homogeneous form via the < p{r)^ > average of the CDA 
method. Given a gaussian distribution, 

^ ^dV|<A(r)n</>(r)p 



(31) 



We define the effective Ts parameter corresponding to 
the CDA density by p = l/(7rf2). 



(32) 



The correlation energy in the CDA, E^^°- for the inho- 
mogeneous distribution, inclusive of spin-polarization ef- 
fects, is now evaluated using fs in any of the well-known 
2D-functionals[3- Note that for typical values of a at 
Ts = 20, the CDA density parameter is ~ 400, while at 
rs = 100, it becomes --7000. Thus we see that the CDA 
replaces the inhomogeneous fluid with sharp Gaussian 
peaks by a uniform, ultra-low- density 2D fluid. In calcu- 
lating Ec{fs) using, say, the formula due to Attaccalite et 
al., a difficulty arises since it is fitted to a maximum of 
40, together with asymptotic forms, while the CDA calls 
for Ts values which are one or two orders bigger. Never- 
theless, we find surprisingly good results (see below). 

At this point we ask if the exchange energy, evalu- 
ated for this ultra-low density fluid, should also be in- 
cluded. We believe that this is indeed the case. The 
fixed-gaussian Wigner-exchange, Ea. l3UI simplv does not 
allow any exchange, and ignores the possibilities of tun- 
neling, ring-exchange etc., that exist in the system. We 
consider that the estimate of exchange obtained from the 
ultra-low density fluid of the CDA accounts for such ex- 
change effects. This point of view is justified post-facto by 
the good agreement of our total energies with the QMC 
total energies. 



C. The LDA exchange and correlation energies 

Another approach to replacing the inhomogeneous 
electron density by a homogeneous fluid-density is the 
local-density approximation (LDA) 22]. Here a uniform 
density corresponding to each local density n{r) is in- 
voked. Thus a local-density parameter fs(r) is defined 
by 



p{r) 



(r) 



(33) 



Hence, knowing the exchange-correlation energy density 
Cxc for a homogeneous system, the exchange-correlation 
energy of the inhomogeneous system is given by 



TABLE III: Results of energy minimization for a hexagonal 
lattice. The Madelung energy Em = — l.lOGlOS/rj, has been 
subtracted out from the total energy. The CDA and LDA 
results are compared with the GFMC calculations of Tanatar 
& Ceperley, T-C The energies are in 10~^ atomic units. 



rs 


20 


30 


40 


50 


60 


CDA 
LDA 
T-C 


0.9247 
0.9404 
0.9167 


0.4824 
0.4899 
0.4983 


0.3059 
0.3102 
0.3234 


0.2156 
0.2185 
0.2313 


0.1625 
0.1645 
0.1758 



Just as in the CDA, the LDA demands the evaluation 
of Ec at densities which are beyond the range of the 
standard fits. Thus LDA needs rs{r) - 300 to 5000 at 
Ts = 20, i.e., a little less extreme than the CDA. Hence, 
some of the shortcomings of the LDA may also be due 
to poorly known correlation energies at the exceptionally 
high values that are required. The LDA can be further 
improved by including gradient corrections. However, we 
have not included them in this study. 



D. Minimization of the total energy Et 

We have now all the energy contributions needed 
to calculate the ground-state energy of a perfect two- 
dimensional (i.e., J-thin) Wigner crystal at a given value 
of the density parameter rs. The energy minimum with 
respect to a is found to be insensitive to the polar- 
ization of the lattice. This in agreement with previ- 
ous studies 0, 0, S El. In Table Hn we give the 
energy correction to the Madelung energy obtained by 
the minimization of Et, using the CDA or the LDA 
for evaluating the exchange-correlation effects, together 
with the results of previous work QMC results by 
Rapisarda and SenatoreQ are very similar to those of 
Tanatar et al., and the agreement is similar. The opti- 
mal a which minimizes the energy is found to be given 
by a = a/rf/^ with a=0.639 for both CDA and LDA ap- 
proaches. Here we note that Nagy's method 18] predicts 
a value a„g = 0.5/rs^^ as a lower bound. A crucial test 
of the accuracy of the CDA and LDA would lie in their 
ability to predict the liquid— )-solid phase transition. This 
is addressed in section IIV CI The total energy can be 
represented by: 



0-2 



as 



,V2+f2+0{rs~'/') r.»l (35) 



E. 



LDA 



= / d'^rexc[rs{r)]p{r) 



(34) 



where ai = —1.106103 is the Madelung constant and 02 
is the zero-point energy of the lattice. We determined 
the coefficient 03 by a least-square fit. The results are 
summarized in Table Hvl together with previous results. 
These results justify our use of the CDA and the LDA 
for evaluating the total energy of quasi-2D Wigner crystal 
phases for which there are no QMC calculations as yet. 
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TABLE IV: Coefficients ai - 03 in Eq.gSlfitting the CD A and 
LDA total energy (for the range rs— 20 to 100) are compared 
with previous work. 





CDA 


LDA 


BMf7] 




TCf5] 


— ai 


1.1061 


1.10610 


1.1060 


1.104715 


1.10610 


02 


0.8142 


0.8142 


0.8142 


0.7947 


0.8142 


as 


0.2456 


0.1194 




0.07338 


0.0254 



IV. INFLUENCE OF THE THICKNESS 

We consider a quasi-2D electron crystal where each 
electron is localized at each lattice site with a gaus- 
sian distribution centered on each site in the {x, 2/)-plane, 
while the z-extension may typically have the form of a 
Fang-Howard density. As before, such z-distributions can 
be replaced by a constant-density form for ease of cal- 
culations. Also, we assume that the 2D layers are in 
HIGFETS, and as such the FH-parameter b is automati- 
cally specified (via Eq.|^)whcn the rs parameter defining 
the 2D-layer density is specified. 

The kinetic energy and the harmonic energy of the 
quasi-2D system are still given by Exia) — Na since 
this is a result of the assumed gaussian form of the wave- 
function. However, the simple Coulomb potential 1 /r has 
changed to F{r)/r where F{r) is the form factor arising 
from the subband distribution. The Wigner-exchange 
energy, i.e., the exchange between two localized electrons 
is now even weaker than in Eq. Hence this type of 
exchange is totally negligible. 



The evaluation of E^c for thick-2D layers using 
CDA and LDA. 



As described in Eq. |51 the z-distribution is mapped 
onto a uniform slab of width a; in HIGFETS this is 
directly related to the parameter in the 2-D plane. 
The inhomogeneous 2-D distribution in the plane is also 
mapped onto a homogeneous distribution via the CDA 
as in Eqs. EHor as in Eq. 01 for the LDA. Both CDA 
and LDA require a knowledge of the Exc{ts, o) for quasi- 
2D uniform systems with a layer width a. Parametrized 
forms of the exchange energy and the correlation energy 
of uniform quasi-2D electron fluids have been given by 
Dharma-wardanafioj . The exchange energy (r-j , a) is 
given in the form £'a;(rs, C)(5(''s, a, C) where Ex{rs,C) is 
the well-known exchange energy of the 5-thin system, 
while Q{rs, a, C) is a form factor. The correlation energy 
of quasi-2D layers in HIGFETS is given in ReL ^ as an 
interpolation involving a form for electron "rods" inter- 
acting via a logarithmic potential (as is the case for small 
Ts), and for 3D-like electrons when rs, and the thickness a 
become large. The Wigner crystallization regime involves 
the latter regime. For details of these parametrizations, 
the reader is referred to ReL .10]. Since the WT involves 
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FIG. 3: (Color online) Comparison of liquid and solid-phase 
energies. {E - EmYs^^'^ where Em = -1.106103/r3 and E 
is the unpolarized or fully-polarized fluid energy, or the solid 
energy E^da, Eua or from QMC. 



small energy differences, we have actually done an ex- 
plicit calculation instead of using the fits. 



B. Minimization of the total energy Et 



As in Sec lIII Dl we minimize the total energy as a func- 
tion of a for a given r^. Here, the energy is more sensitive 
to the spin-polarization C, than in the perfect crystal even 
if the difference is very small. The unpolarized crystal is 
more stable than the polarized one. So we present results 
for the unpolarized system. The values of a which mini- 
mizes the energy can also be fitted by the same form as 
in Sec llHDl We obtain 



Otcda 



0.619 

'ZJJ2 



and aida = 



0.627 

J, 3/2 



(36) 



We have also fitted the total energy. Here the Madelung 

energy is the Ecou given in Eq. H25|) and we use the usual 

• • 3/2 , 

expansion m inverse rs etc. 



jpcda 

H/rp 



^ , . 0.68597 0.321652 

Ecou {rs) + ^ + ^ (37) 



0.708977 0.357242 



3/2 



(38) 



We remark that the total energy has a minimum as a 
function of r^. This minimum is located around ^ 26 
and its value is ^ —0.011 a.u. A comparison of liquid 
phase and Wigner crystal in HIGFETS is given in Ta- 
ble These total energies include the 2a/3rs correction 
arsing from the interaction of the quasi-2D layer with 
the unifrom background, as discussed in subsection III Bl 
Since this depends on the layer thickness a, this correc- 
tion does not occur in the ideal 2D system. 
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TABLE V: Results of energy minimization for a hexagonal 
lattice and comparison with the unpolarized liquid phase en- 
ergy El- The energies are measured in 10""^ atomic units. 



rs 


15 


20 


30 


50 


Ecda 


-6.7255 


-10.1581 


-10.8576 


-9.1112 


Elda 


-6.5036 


-9.8169 


-10.7782 


-9.0306 


El 


-7.1324 


-10.0249 


-10.5995 


-8.8939 



Wigner transition occurs at = 15.5 for the CDA, and 
Ts — 20.5 for the LDA, i.e, before the spin-phase tran- 
sition (marked SPT in the figure) of the liquid phase. 
Since the (5-thin 2D layer is expected to have a WT near 
Ts ^ 35, the thickness effect has brought the WT to 
smaller Vg values. It should be noted that if correla- 
tion effects are neglected, the WT occurs at very low Vg. 
Hence the shift of the WT to low Vg is a consequence of 
the reduced correlations in the quasi-2D system. 




15 20 25 

r 

s 



FIG. 4: (Color online) {E - EmYJ^'^ where Em = 
— 1.106103/rs and E is the unpolarized or fully-polarized fluid 
energy, or the solid energy Ecda or Eida- The spin-phase tran- 
sition in the fluid is labeled SPT. The onset of the Wigner 
crystal in CDA and LDA are indicated by arrows. 



V. CONCLUSION 

We have investigated the crystallization of electrons 
under their own interaction in infinitely-thin 2D electron 
layers, known as Wigner crystallization, as well as in 2D 
layers with a finite width. The case of infinitely-thin 2D 
electron layers has been extensively studied in the past, 
and provides a bench mark to test our methods for re- 
placing inhomogeneous electron densities by equivalent 
uniform-density models. Detailed Hartree-Fock calcula- 
tions with large plane-wave basis sets are shown to be 
closely equivalent to the single-gaussian harmonic lat- 
tice calculations. We show that the constant-density ap- 
proximation, CDA, as well as the local-density approx- 
imation, LDA, successfully recover the total energy as 
well as the correlation energies of the infinitely-thin 2D 
electron crystal. In particular, these models predict a 
liquid— >solid phase-transition in the range 30 < rg < 40, 
in good agreement with Quantum Monte Carlo simu- 
lations. When these methods are applied to quasi-2D 
layers with the thickness as in HIGFETS, the weak- 
ened Coulomb correlations move the Wigner transition 
towards high densities. The LDA and the CDA predict 
a Wigner transition near rs ~15 to 21. 



C. Phase transition Liquids Wigner crystaL 

According to quantum Monte Carlo simulations, the 
phase transition between a (5-thin 2D electron liquid and 
a 2D electron Wigner crystal occurs around r<; = 35 ± 5. 
With our methods we find a transition for — 32 using 
the CDA and = 38 using the LDA. 

In Figure 131 we show the phase diagram of the system 
(in order to have a clear display we present {E—EM)rs^^'^ 
where Em — — 1.106103/rs the Madelung energy). The 
fluid phase energy is calculated using the fit given byQ. 
Unlike in the i5-thin 2D system, the total energy contains 
the term A^e = 2a/3r^ arising from the interaction with 
the unifrom background, (see subsection llllj|l . This has 
been removed from both the liquid and the solid phase 
energies as this improves the clarity of the display. 

These results tend to show that both LDA and CDA 
provide an adequate evaluation of E^c for electron solids, 
especially when we recall that the Ec at values (200- 
8000), way outside the fit range (r^ < 40), are needed in 
the CDA and the LDA evaluations. Figure 01 displays the 
phase transitions in the quasi-2D HIGFET system. The 



APPENDIX: EVALUATION OF THE 
ELECTROSTATIC ENERGY AND LATTICE 
SUMS. 

The expression for the electron-electron Coulomb en- 
ergy. Eg. 1241 can be rewritten, using the Jacobi function 
technique as. 



Eeeiq) 



2.^/11 ^ Jo 



-v\pi\ -»q-pj 



yaa 



where J^.o is the Kronecker symbol, because when pj = 0, 
we must have (m, n) ^ (m', n'). 

The advantage of introducing the Jacobi function 9{z, X) 
is that it converges well for large X and we are also able to 
obtain convenient well-convergent formulas for the small- 
X region by applying the transformation fEa. l20() ) from 
which the Coulomb singular part at q — > can be rigor- 
ously extracted. Thus, Eee{q) obtained in Ea. (|A.l|l can 
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be separated into a large y part and a small y part given where 

by 



(A.l) 



and 



, Jo 



-y\pj\ — iq-ft 



2aia2 — JO 

ya 



N 



E / dyy-^/^me-i^ ^ 



n 



n ' 

\a—x,y 



r d,y-'"me- 

2aia2 ^ Jo 



2aia2 



dyy-^/'^,f(y)e-^ 



Vo 



n 



\a.—x,y 



^{2p^y + ^q^)■ 



27r 



2tt 



{2pfy + iq^) 



2pfy + iqa n 



2aay ' ya^^ 



5.U 



1 + 1 



\a—x.y 



- 1 



20? 7o 



2aay yao 
dyy-'/'f{y)+E':r{q) 



(A.2) 



(A.3) 



where n/ is the number of sublattices (0102/71; = Trrg^). 
The results of the calculation are independent of the 
value of yo > 0; nevertheless, we choose it such that 
the sums 9 converge fast, and we have 



E. 



>-"^iq) = !^l dyy-'/'f{y)e-ij (A.4) 
zai02 



In order to complete the calculation of Ecom we need to 
discuss the form of In their article (l&j . Fujiki and 

Geldart use the Fang-Howard density ?7//t(^) = ^z^e~'s 
(Figure^. As already discussed we replace the FH dis- 
tribution by the equivalent CDA, i.e., we use a constant 
density slab of width a = 16&/3. With this homogeneous 
form of density 

/•a pa 

f{y) = f dz dz'e-y^'-'">^ 

Jo Jo 
^ _2 (e-°'" + oV^erf(oVy)-l) 

y 

We see here an advantage of the constant density map- 
ping to density 77, the analytic expression of f{y) being 
quite simple. 



In Ea. l|18|l . we replace r]{z') by its expression 



Eed{q) 



-{l-e-'"i) 



i?ed(g-0) --—(-- --6, + 0(g)) (A.6) 



qrs aq 
N fl a 
Ts^ \q 2 

We recall that bd is positive or zero, and gives the loca- 
tion of the donor layer at z — — 6^- Now, in Ea. (jA.4p . we 
use the definition of f{y)- 



E!:r{q) = 



N 2 



1 



qrs'' qa'' \ q 



^^ee°'"(9->0) = ^Q-^+O(,) 



(A.7) 



Now, we will use Eq.lISl), jSS, JSSl), 

in Ea. (|12|l to calculate the Coulomb energy for different 
types of lattices and for different thicknesses. We can 
see that the expression of Ecou is dependent on the pa- 
rameter bd- Since this is a constant contribution, we set 
bd = and focus on the part which depends only on the 
geometry of the lattice and on its thickness. 
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